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ABSTRACT 

New  algorithms  are  described  for  Bayesian  estimation  of  parameters 
in  nonlinear  models  of  multiple-response  systems.  Modal  and  interval 
estimates  are  provided  for  the  parameter  vector  9  of  the  predictor 
model,  and  for  the  variance-covariance  matrix  c  of  a  Normal  error 
distribution.  Allowance  is  made  for  gaps  (missing  values  of  responses) , 
such  as  commonly  occur  in  practice.  Two  chemical  examples  are  analyzed. 
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SIGNIFICANCE  AND  EXPLANATION 


Some  new  algorithms  are  presented  for  fitting  mathematical 
models  to  multiple-response  experiments.  These  algorithms  give 
estimates  of  the  parameters  in  a  user-defined  predictor  model, 
and  also  estimate  the  parameters  of  a  Gaussian  model  of  the 
observational  error  distribution.  The  development  is  based  on 
Bayes*  theorem,  and  provides  a  natural  extension  of  known  least- 


squares  estimation  methods.  Allowance  is  made  for  missing  values 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


NEW  ALGORITHMS  FOR  NONLINEAR  LEAST  SQUARES 
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New  algorithms  are  described  for  Bayesian  estimation  of  parameters  in 
nonlinear  models  of  multiple-response  systems.  Modal  and  interval  estimates 
are  provided  for  the  parameter  vector  0  of  the  predictor  model,  and  for  the 
variance-covariance  matrix  a  of  a  Normal  error  distribution.  Allowance  is 
made  for  gaps  (missing  values  of  responses),  such  as  commonly  occur  in 
practice.  Two  chemical  examples  are  analyzed. 

INTRODUCTION 

Realistic  models  of  multivariate  phenomena  often  relate  several  predicted 
responses  to  a  common  set  of  parameters.  Multiresponse  experiments  are  re¬ 
quired  to  establish  such  models,  but  frequently  yield  irregular  data  which 
are  difficult  to  analyze  by  classical  methods. 

Bayes '  theorem  is  a  good  starting  point  for  parameter  estimation  in  these 
situations.  The  multivariate  error  distribution  can  be  estimated  concurrently, 
whereas  it  has  to  be  prescribed  when  least-squares  methods  are  used.  Thus, 
the  Bayesian  approach  allows  more  objective  parameter  estimates,  if  sufficient 
data  are  provided.  An  excellent  general  account  of  this  approach  is  given  by 
Box  and  Tiao  (1973) . 

Bayesian  inference  deals  with  a  data  array  {y  =  y,  a  model  for  E(y) 
with  parameter  vector  0  ,  and  an  error  distribution  model.  If  a  Normal  error 
model  is  used,  with  variance-covariance  matrix  <j  ,  the  unknown  elements  of 

c  will  appear  as  additional  parameters.  The  full  set  of  parameters  can  be 
estimated  optimally  by  maximizing  the  posterior  density  p(6,o|y)j  confidence 

regions  can  also  be  calculated  from  this  function. 
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In  certain  cases,  the  posterior  density  can  be  integrated  analytically 
to  obtain  the  marginal  density  p(Ojy).  Box  and  Draper  (1965)  accomplished 


this  for  multivariate  Normal  error  distributions  and  rectangular  data  structure 
(Table  la).  For  block-rectangular  structures  (Table  lb),  p(9|y)  is  the  prod¬ 
uct  of  the  Box-Draper  densities  for  the  individual  rectangles.  More  compli¬ 
cated  data  structures  often  occur,  however,  such  as  that  in  Table  lc,  for 
which  p(0 |y)  cannot  be  expressed  in  closed  form.  Therefore,  in  this  paper 
we  use  the  full  posterior  density  p(6,oJy),  which  lias  a  closed  form  for  any 
finite  data  structure. 

Inspection  of  the  parameter  estimates  and  residuals  often  suggests 
alternatives  to  the  postulated  model.  Therefore,  parameter  estimation  should 
not  be  viewed  as  an  end  in  itself,  but  should  be  followed  by  critical  examina¬ 
tion  of  the  model  and  investigation  of  any  promising  alternatives.  Interesting 
predictions  or  unresolved  differences  between  models  will  naturally  lead  to 
further  experiments . 

Table  1.  Examples  of  Data  Structures  with  m  =  4  and  n  =  8 


u 

1 

2 

3 

4 

5 

6 

7 

8 


la.  Rectangular 

yul  yu2  yu3  yu4 

+  +  +  + 

+  +  +  + 

+  +  +  + 

+  +  +  + 

+  +  +  + 

+  +  +  + 

+  +  +  + 

+  +  +  + 


lb.  Block-rectangular 

yul  yu2  yu3  yu4 

+  + 

+  + 

+  + 


+  + 

+  + 

+  + 

+  + 

+  + 


lc.  Irregular 

*  ul  yu2  yu3  \A 

+  +  + 

+  +  +  + 

+  + 

+ 

+  + 

+ 

+  +  + 

+  + 


PROBLEM  FORMULATION 


Consider  a  set  of  independent  experiments,  u  =  l,...,n,  in  which  a 
table  {y  of  observed  responses  have  been  obtained  at  known  settings 
(x^l  of  the  independent  variables.  There  are  m  linearly  independent 
kinds  of  observations;  thus  the  index  i  ranges  from  1  to  m  ,  but  in 
each  experiment  some  values  may  be  missing  as  in  Tables  lb  and  lc. 

The  observations  in  the  uth  experiment  are  regarded  as  a  sample  from 
a  population  of  the  form 

(1) 

The  functions  f^(xu,0)  are  models  for  the  expected  responses  E(y^.  1  (*)  • 
The  residuals  e  ^  in  the  uth  experiment  are  treated  as  a  random  sample 
from  an  m-variate  Normal  distribution;  this  gives  the  probability  density 
(Wilks,  1962) 

-m  /2  , 

n  ,  .-1  /9  t  -1 

(2) 


y  .  =  f . (x  ,6)  +  e  . 
ux  l  ~u  -  ui 


p(€u|o)  =  (2ir)  u  |o  |  1/2  exp(-J  cT*  cj. 


-u  -  u  ~u 

Here  e  is  the  column  vector  of  error  variables  e  ,,...,  e  with  dummy 
u  ul  um 

zeroes  inserted  where  observations  are  missing.  Correspondingly,  o  is 

obtained  from  the  full  variance-covariance  matrix,  o  =  (o  }  ,  by  sub- 

13 

stituting  dummy  elements  6.  .  whenever  observation  y  .  or  v  is 

13  ui  Juj 

missing.  Here  5„  is  unity  when  i=j,  and  zero  otherwise. 

The  joint  error  density  model  for  the  set  of  n  experiments  follows 

directly  from  Equation  (2) ; 

n  -m  /2 

p(c|e,o)  =  n  (2tt)  u  |oJ'1/2  exp (-£  ET  o"1  e  ).  (3) 

-  ~  ~  u=i  ~u  'u  ~u  -u 

Insertion  of  Equation  (1)  gives  the  corresponding  density  in  observation 
space : 


n  -m  /2 

P<y|e,o)-[n  (2;)  u  \a\  /2) 

u=l  “u 


(4) 


m  m 


exp{-J  l  l  l  o*3  [y  .  -  f  .  (0)]ty 

u=l  i=l  j=l  u  U1  U1  "  u3 


V2)1}. 
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Here  the  functions  fu^(0)  stand  for  f^(xu,0)  evaluated  at  the  known 
settings  x^  of  the  independent  variables.  The  a*3  are  the  elements  of 
the  precision  matrices  o^1.  The  right-hand  term  may  also  be  regarded, 
by  Bayes'  theorem,  as  the  likelihood  function  for  0  and  o  when  eval¬ 
uated  with  given  observations  y  . 

The  usual  factorization  of  the  prior  density  p(0,o)  is  assumed, 
p(6,o)  =  p ( 8 )  p(g)  (5) 

and  a  locally  uniform  density  p(0)  is  assumed  in  the  region  of  appreci¬ 
able  likelihood.  The  latter  assumption  requires  some  care  in  the  parame- 
trization  of  the  model.  The  prior  density  of  o  is  taken  from  Box  and 
Draper  (1965) : 

P(o)  •*  |o|-(m+1,/2  •  (6) 

Bayes'  theorem  then  gives  the  posterior  density 


p(0,a|y)  =  p(0,g)  p(y|0,n) 

=  c|cr(m+1,/2 1  n  ig  r1/2i 

_ i  ^ 


(7) 


n  m  m 

•  exP{-i  l  l  T  a13ly  .  -  f  . (0)]  [y  .  -  f  .  (0)]} 

.  \  ,L.  u  ui  ui  -  J  un  UT  - 

U=1  1=1  ]=1  J  J 

in  which  c  is  a  proportionality  constant.  All  that  the  data  reveal  about 

the  parameters  0  and  o  is  contained  in  this  density  function. 

Point  estimates  of  0  and  a  are  obtainable  by  maximizing  the  posterior 

density  just  described,  or  by  minimizing  the  function 

S(ij>)  =  S(0,a)  =  -2  In  p(0,a|y)  +  2  In  c 

n~  ~  ~ 

=  (m+1)  In  J  a  |  +  £  In  |a  |  (8) 

u=l  ~u 

n  m  m  , 

+  I  l  l  of,3  (y  -  f  .  (9)  1  [y  .  -  f  .  (9)) 

u«l  i-1  j=l  u  ui  U1  -  U3  u3  > 

over  the  permitted  region  of  0  and  a  .  Here  t|>  is  a  column  array  of 
the  model  parameters  0^,  ...,0^  and  the  independent  elements  of  o  . 


> 


I 


I 


\ 
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The  latter  are  taken  from  the  lower  triangle  of  a  in  row  order,  i.e. 
i|<  v  =  a.  .  with  i  >  j  and  k  =  j  +  i(i-l)/2.  Thus,  the  total  number 

p+K  — 

of  parameters  is  q  =  p  +  m(m+l)/2. 

If  the  matrix  a  were  believed  to  be  known,  i.e.,  if  a  sharply  focussed 
prior  density  p(g)  were  assumed,  then  S ($)  would  reduce  to  S(0)  and  we 
would  have  a  least~squares  estimation  problem  with  just  p  parameters.  In 
practice,  one  seldom  knows  o  accurately;  hence,  the  full  Bayesian  solution 


is  recommended. 


PARAMETER  ESTIMATION  ALGORITHMS 

Several  algorithms  are  described  here  for  obtaining  summary  informa¬ 
tion  from  Equation  (8) .  These  algorithms  are  part  of  a  Fortran  IV  package 
available  from  the  authors. 

1.  Counting  Algorithm 

Before  analyzing  S  we  count  Equations  (1)  to  see  which  parameters 

can  plausibly  be  estimated  from  the  data.  We  first  try  to  match  each 

parameter  o,  .  in  b  with  an  observation  pair  (y  of  a  replicate 

kg  I  uk  ug 

experiment  (i.e.,  an  experiment  which  has  the  same  expected  response  values 

as  a  prior  experiment  in  the  data  set) .  If  this  process  cannot  be  completed 

for  a  given  k  ,  we  then  try  to  match  each  remaining  error  parameter  o.  .  , 

kg 

and  each  model  parameter  0  in  the  function  pairs  If  ,  (0),  f  .(0)1,  with 
r  r  ^  uk  -  uj  - 

a  non-replicate  observation  pair  ^y^y  )  •  Finally,  any  remaining  model 
parameters  0^  are  matched  with  remaining  non-replicate  observations.  If 
the  matching  can  be  completed  for  all  elements  of  0  ,  we  proceed  with  the 
estimation.  Otherwise,  the  full  set  of  parameters  cannot  be  estimated  from 
the  data. 

The  counting  algorithm  is  a  logical  Gaussian  elimination .  This  test 
is  a  useful  diagnostic,  but  is  not  infallible,  since  the  actual  rank  of  the 
estimation  equations  depends  on  the  numerical  values  of  x,  y,  and  0 . 

2.  Minimization  Algorithm 

A  modified  Newton  method  is  used  to  find  a  minimum  of  S(iJ;).  Let  tJjQ 
be  the  value  of  \p  at  the  start  of  an  iteration.  A  correction  vector 
(^1  ~  V  comPu*'c<^  by  minimizing  the  local  quadratic  expansion  (see 

Appendix  A  for  derivative  expressions) 


2 

tT  3  S 


s(!)  =  s{y  +  0  (!  •  y  +  J  -  y  o(*  -  y 
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over  a  user-specified  rectangular  region  around  .  The  region  is  chosen 

small  enough  to  ensure  that  S(^)  is  a  good  approximation  to  the  function 
S(iJ')  of  Equation  (9).  A  search  is  then  made  for  a  minimum  of  S  in  the 
interval  of  positive  definite  o  on  the  line  from  ^  through  4^  ;  this 
gives  the  starting  point  for  the  next  iteration.  The  calculation  continues 
until  two  successive  line-minima  agree  within  confidence  intervals  calculated 
from  Equation  (14)  for  each  parameter. 


3.  Response-Independence  Test 

Box  and  co-workers  (1973)  have  pointed  out  the  need  to  test  the  responses 

for  linear  independence.  Preferably,  one  should  perform  this  test  on  the 

residuals  fy  .  -  f  (0)],  which  miqht  become  linearly  dependent  in  certain 
ui  ui 

regions  of  0  .  In  the  present  procedure,  such  linear  dependence  is  readily 

detected  during  the  inversion  of  o  at  the  start  of  each  iteration.  The 

calculation  can  continue  if  all  pivot  elements  (Stewart,  1973)  found  in  this 

inversion  are  greater  than  a  specified  fraction,  say  0.1,  of  the  corresponding 

elements  o. .  . 

n 

4.  Confidence  Regions 

Equation  (8)  gives  the  simple  form 

p(4>|y)  a  expt-i  S(*)l  (10) 

for  the  posterior  density  function,  or  "confidence  density".  Use  of  liquation 
(9)  gives  the  approximation 

P ( 4* | y )  “  exp[-i(4*  -  4,)T  A(4>  -  (11) 

valid  in  the  neighborhood  of  the  minimum  point  4'  •  Here  A  is  the  qxq 
matrix  (positive  definite  since  S  is  at  a  minimum)  with  elements 


(12) 
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computed  as  described  in  the  Appendix.  Thus,  near  the  optimum,  the  param- 

—  1 

eters  are  Normally  distributed  with  variance-covariance  matrix  A  .  If 
Equation  (11)  is  used  as  an  approximation  for  all  values  of  'y  ,  then  the 
confidence  intervals  for  Normal  distributions  can  be  applied.  For  example, 
the  ellipsoidal  region 


roughly  approximates  the  100(1  -  a)  percent  highest-posterior-density  region, 
or  joint  confidence  region,  for  ip  based  on  the  given  data.  The  intervals 

(j<^k  “  k\/  h  Akk  )  <  erfc-1  (a)  (14) 

roughly  approximate  the  100(1  -  a)  percent  confidence  intervals  for  the  indi¬ 
vidual  parameters.  For  symmetric  95  percent  confidence  intervals  (a  =  0.05)  , 
erfc  1(a)  has  the  value  1.96. 

Equation  (14)  is  more  reliable  than  (13),  since  the  integration  used  to 
obtain  it  is  less  affected  by  the  tails  of  the  posterior  density  function. 

More  accurate  intervals  can  be  obtained,  but  with  greater  effort,  by  numerical 
integration  of  Equation  (7)  or  (10). 
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RESULTS  FOR  RECTANGULAR  DATA  STRUCTURES 


(21) 


Equation  (20)  gives,  at  the  stationary  point, 

v  (6) 

*  rs  ~ 

o  =  - : —  . 

rs  m+n+1 

Hence, 

orS  =  (m+n+1)  vrS(0)  .  (22) 

Insertion  of  Equation  (22)  into  (19)  gives  Equation  (18)  at  the  stationary 
point  of  p(G,o|y).  Hence,  for  rectangular  data  structures,  the  same  values 
of  0  and  o  are  obtained  whether  one  maximizes  p(9,o|y)  or  p(0|y). 

Of  course,  the  marginal  confidence  regions  for  0  can  be  estimated  more 
directly  in  the  latter  case.  The  normal  equations  based  on  p(0|y),  given 
by  Stewart  and  Sorensen  (1976) ,  are  convenient  for  this  purpose. 

The  covariance  estimates  in  Equation  (21)  are  maximum-density  values, 
and  thus  differ  from  the  expectation  values  Efo^Jy)  unless  n-m-p  is 
very  large.  If  expectation  estimates  of  the  are  desired,  one  can  com¬ 

pute  them  as  the  corresponding  moments  of  the  normalized  posterior  density 

p(e,o|y) . 

EXAMPLE  1.  Kinetics  of  a  Three-Component  S’,  stem 

Consider  the  chemical  conversion  of  initially  pure  species  1  to  species 
2  and  3  in  a  batch  isothermal  reactor.  Simulated  data  for  the  system  are 
given  in  Table  1,  reproduced  from  Box  and  Draper  (1965);  here  y  is  the 
yield  of  species  i  in  experiment  u  .  The  system  is  modelled  by  the  differ¬ 
ential  equations 


10- 


fx  =  expf-^  t) 


f2  =  (expt-^  t)  -  exp(-k2  t))kj/(k2  -  k^ 
f3  =  1  "  fl  -  f2 

under  the  indicated  initial  conditions.  As  noted  by  Box  and  Draper,  it  is 
natural  to  regard  the  parameters  8^  =  In  k^  as  uniformly  distributed  a^  priori . 

There  are  three  responses  y  ^  per  experiment.  Only  two  would  be  linearly 
independent  if  the  yields  were  mass-balanced  (i.e.,  if  the  yields  in  each  row 
added  up  to  unity).  The  data  in  Table  2  are  clearly  not  mass-balanced,  so  we  use 
all  three  columns  of  responses. 

The  replicates  in  Table  2  allow  preliminary  estimation  of  the  parameters 
a.  by  the  relation 

IT 

J  n 

R 


1 

’ij  2n  ^  ^yri 
J  R  r=l 


-  y ' . )  (y  . 
n  ri  rj 


y’  .)  • 
*1 


Here  y  ^  and  y^  are  the  observations  of  response  i  in  the  first  and 
second  tests  of  replicate  pair  r  ,  and  n  is  the  number  of  such  pairs. 
This  procedure  gives 

0.00102  -0.00128  0.00025 

-0.00128  0.00351  0.00024 

0.00025  0.00024  0.00101 


(Si.)  - 


as  a  preliminary  expectation  estimate  of  o  .  This  is  a  well-conditioned 
matrix,  so  our  choice  m  =  3  was  correct. 

The  parameter  vector  i p  for  the  present  example  consists  of  0^,  0^, 
and  the  six  elements  on  and  below  the  diagonal  of  o  .  To  test  the  conver¬ 
gence  of  the  estimation  from  a  poor  initial  guess,  the  calculation  was  started 
from  the  initial  value  shown  in  Table  3.  Convergence  was  obtained  in  eight 
iterations,  to  the  point  estimates  and  95  percent  confidence  intervals  given 
there. 
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Table  2.  Data  for  Example  1,  from  Box  and  Draper  (1973) 


tu 

»“H 

>% 

yu2 

yu3 

0.5 

0.959 

0.025 

0.028 

0.5 

0.914 

0.061 

0.000 

1. 

0.855 

0.152 

0.068 

1. 

0.785 

0.197 

0.096 

2. 

0.628 

0.130 

0.090 

2. 

0.617 

0.249 

0.118 

4. 

0.480 

0.184 

0.374 

4. 

0.423 

0.298 

0.358 

8. 

0.166 

0.147 

0.651 

8. 

0.205 

0.050 

0.684 

16. 

0.034 

0.000 

0.899 

16. 

■  0.054 

0.047 

0.991 

Table  3. 

Parameter 

Values 

for  Example  1 

Initial 

Solution  1 

Solution  2 

Solution  3 

Parameter 

Value 

Eqs.  (8.14)* 

Eqs.  (8,14)* 

Eqs. 

(18,21)* 

81 

-2.3026 

-1.5723+0.0567 

-1.5723+0.0558 

-1.5723+0.0800 

°2 

0. 

-0.7023+0.1374 

-0.7023+0.1346 

-0.7023+0.1931 

-3 

(0.76+0.53) 

-3 

-3 

°11 

0.01 

(0.76+0.52) 

10 

10 

0.76 

10 

-3 

-(0.50+0.63) 

-3 

-3 

°21 

0. 

-(0.50+0.63) 

10 

10 

-0.50 

10 

-3 

(1.86+1.29) 

,  -3 

-3 

°22 

0.01 

(1.86+1.28) 

10 

10 

1.86 

10 

-3 

(0 . 32±0 .41) 

-3 

-3 

°31 

0. 

(0.32+0.41) 

10 

10 

0.32 

10 

-3 

(0.40+0.62) 

-3 

-3 

°32 

0. 

(0.40+0.62) 

10 

10 

0.40 

10 

°33 

0.01 

(0.77+0.54) 

io“3 

(0.77±0.54) 

10-3 

0.77 

io-3 

* 

All  intervals  are  95%  highest  posterior  density  regions.  In  solution  3,  the 
intervals  are  computed  from  the  normal  equations  with  ’’residual  mean  square" 

A 

|v(0)|/(n-2)  and  n-2  *  10  residual  degrees  of  freedom.  In  Solution  1,  the 
second-derivative  terms  of  Equation  (AlO)  are  included. 
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A  second  calculation  was  made  with  the  same  initial  values,  but  with 


second-order  O-derivatives  neglected.  Convergence  was  obtained  to  the  same 
point  estimates  in  nine  iterations.  The  confidence  intervals  differed 
slightly,  as  shown  in  Table  3. 

A  third  calculation  was  made  by  minimizing  the  determinant  ]v(6)  |.  Box 

and  Draper  (1965)  did  this  by  a  search  procedure;  we  used  the  modified  Newton 

algorithm  of  Stewart  and  Sorensen  (1976) ,  but  neglected  the  second-order  0- 

derivatives  of  the  functions  f  , (0)  .  Convergence  was  obtained  in  seven 

ui  _ 

iterations,  tc  the  same  point  estimates  9^  .  The  point  estimates  for  the 

,  computed  from  Equation  (21) ,  also  agreed  exactly  with  the  two  preceding 
solutions.  The  one-parameter  confidence  intervals  (computed  in  this  case  only 
for  0^  and  02)  are  wider  than  before,  and  are  considered  more  accurate 
since  in  this  case  the  o_  have  been  integrated  out  exactly  (Box  and  Draper, 
1965). 

EXAMPLE  2.  Kinetics  of  a  Five-Component  System 

Fuguitt  and  Hawkins  (1945,  1947)  did  extensive  experiments  on  the  liquid- 
phase  thermal  reactions  of  c-pinene  and  its  decomposition  products.  The 
following  products, in  order  of  boiling  point,  were  identified. 


A. 

a-Pinene 

C10H16 

B. 

a  -  and  8-Pyronene 

C10H16 

C. 

Dipentene 

C10H16 

D. 

alio -Oc imen e 

C10H16 

E. 

Dimer 

C20H32 

The  reaction  conditions  and  yields  are  reported  in  Table  3. 

We  have  normalized  the  yields  to  obtain  exact  mass  balances;  this  makes 
the  yields  linearly  dependent,  and  accordingly  we  have  omitted  species  D  . 
The  remaining  species  are  grouped  as  cumulative  distillation  fractions; 


I 


A,  A+B,  A+B+C,  and  E.  Each  of  these  responses  represents  essentially  the  total 
mass  fraction  distilling  above  or  below  a  particular  temperature.  The  yield? 
of  B  originally  reported  in  tests  1-15  have  been  deleted,  since  they  were 
interpolated  values  rather  than  observations  (Fuguitt  and  Hawkins,  1947;  Box 
and  co-workers,  1973). 

There  are  numerous  gaps  in  the  data.  ct-Pinene  (A)  was  reported  in  experi¬ 
ments  1-16,  but  was  considered  negligible  in  the  remaining  experiments,  pyronener. 
(B)  were  reported  only  in  experiments  16-31;  they  proved  difficult  to  isolate 
except  at  small  concentrations  of  a-pinene.  Only  the  dimer  fraction  (E)  war 
reported  in  the  experiments  with  allo-ocimenc  (D)  or  dimer  (E)  as  feed.  The 
simplified  reaction  scheme  proposed  by  Fuguitt  and  Hawkins  (1947)  implies  that 
o-pinene  (A)  and  dipentene  (C)  would  not  be  formed  in  the  latter  experiments, 
but  that  the  other  three  species  would  be  present. 

The  first  eight  experiments  were  used  for  parameter  estimation  according 
to  Equation  (17)  with  m  =  3  by  Box  and  co-workers  (1973) ,  and  by  the  present  , 
authors  (1976) .  The  full  41  experiments  could  not  be  so  analyzed  because  of 
their  irregular  structure;  therefore,  only  rough  estimates  were  obtainable  for 
several  of  the  reaction  parameters.  With  Equation  (8),  on  the  other  hand,  all 
41  experiments  can  be  analyzed. 

We  postulate  the  following  reaction  scheme, 


with  the  following  differential  equations  for  the  concentrations: 


=  -(kl  +  V  *A 


!!a 

dt 

dt 

■ar  -  ki 


2k5  *A 


=  -k_3  «B  +  k3  ^ 


dij> 

"dt 

dp 


_ D 

dt 


k2  *A  +  k-3  “  k3  "  2k4  +  2k-4  *E 


E  2  2 

■  -  -  =  kr  4  +  k„  6  -  k  ,  4 
dt  5  A  4  D  -4  E 

Here  we  have  assumed  equal  densities  for  the  reaction  mixture  and  all  species. 
The  are  molar  concentrations  relative  to  the  molar  density  of  pure  liquid 

a-pinene  at  the  reaction  temperature.  The  resulting  initial  p.  values  for  the 
pure  reactants  are;  1.0  for  a-pinene,  1.0  for  allo-ocimene,  and  0.5  for 
dimer.  The  rate  coefficients  are  represented  as  Arrhenius  functions. 


In 

<V  -*  ei  - 

(1/T  - 

ei+5 

H* 

II 

H* 

In 

(k3/k-3}  " 

-011/TB 

+  (1/T  -  1/TB) 

°13 

In 

<Vk-4J  = 

-9i2/tb 

+  (1/T  -  1/TB) 

614 

with  k^  values  in  min  ,  T  in  Kelvins,  and  a  base  temperature  Tb  of 


478.5  K. 


The  data  and  parameters  were  paired  to  check  the  feasibility  of  the 
estimation.  This  indicated  a  sufficient  amount  of  data  for  estimation  of  all 
parameters  except  o  ^  .  However,  the  replicate  comparisons  (u  =  18-19,20-21, 
22-23,24-25)  involving  y^2  all  give  duplication  of  y  ;  furthermore  each  of 
these  comparisons  gives  a  duplication  of  either  y^2  or  y  With  these 
results,  we  find  that  neither  a^2  nor  a^2  can  be  estimated,*  indeed,  an 
attempt  to  estimate  them  was  terminated  by  the  linear  independence  test 
described  above.  Thereafter,  021'  °32'  and  °42  were  a*1  zero,  and 

the  remaining  parameters  were  estimated  by  minimization  of  S  . 
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Initial  values  of  the  9-parameters  were  chosen  from  the  results  of 
Fuguitt  and  Hawkins  (1945,  1947),  Box  and  co-workers  (1973),  and  the  present 
authors  (1976).  Initial  variance  estimates  were  calculated  from  repli¬ 

cate  data  available  in  Table  4,  and  zeros  were  inserted  initially  as  covariar. 

The  model  was  integrated,  for  each  experiment,  by  the  method  of  Guerti:. 
et  al  (1977)  with  6  mesh  points.  The  coefficients  in  Equation  (9)  were  com¬ 
puted  as  described  in  the  Appendix,  with  first-order  sensitivities  3$  V  ■  ^ 
computed  by  the  method  of  Stewart  and  Sorensen  (1976) . 

A  first  minimization,  with  reaction  5  omitted,  converged  within  20  iter.: 
tions.  This  gave  S  =  41.06  with  parameter  estimates  as  shown  in  Table  5. 
The  confidence  intervals  show  the  6's  to  be  estimated  quite  precisely.  Tlu 
are  estimated  less  precisely,  as  anticipated  from  the  limited  number  of 
data  on  several  combinations  of  responses.  The  deviations  of  the  data  from 
the  fitted  model  are  shown  in  Table  6. 

A  second  minimization  of  S  was  done  with  the  full  5-reaction  model . 
This  calculation  converged  to  a  very  flat  minimum  at  S  =  34.09,  with  param¬ 
eter  estimates  as  shown  in  Table  5.  The  deviations  of  the  data  from  this 
fitted  model  are  also  shown  in  Table  6. 

The  5-reaction  model  is  better  able  to  describe  the  polymer  yields  from 
a-pinene  at  short  times,  as  can  be  seen  in  Table  6.  We  can  also  test  the 
significance  of  the  added  parameter  6^  by  use  of  the  confidence  intervals. 
Table  5  gives  0,.  =  -11.945  ±  0.698,  based  on  Equation  (14);  this  implies  the 
limits  (1  ±  0.698)  exp(-11.945)  for  k^  with  the  alternate  prior  p(k^)  =  c 
Hence,  the  95%  confidence  interval  for  k^  does  not  include  zero. 

On  the  other  hand.  Equations  (9)  and  (13)  give  the  following  approximate 
expression  for  the  95%  confidence  region  of  the  20  fitted  parameters  of  the 
5-reaction  model; 


■  .tV.: 
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Table  4.  Data  for  Example  2,  from  Puquitt  and  Hawkins  (1945,1947) 


Expt. 

u 

Feed 

• 

O 

e-T 

t  -min 
u 

** 

w 

u 

Normalized  yields,  weiqht 

yul  yu2  yu3 

(A)  (A+B)  (A+B+C) 

percent 

yu4 

(E) 

1 

A 

189.5 

1230. 

1 

88.3 

*** 

96.2 

2.2 

2* 

A 

189.5 

1230. 

1 

88.2 

*** 

95.7 

1.3 

3 

A 

189.5 

3060. 

2 

76.4 

*** 

92.7 

2.8 

4 

A 

189.5 

4920. 

2 

64.8 

*** 

88.9 

5.8 

5 

A 

189.5 

7800. 

2 

50.3 

*** 

84.7 

9.3 

6 

A 

189.5 

10680. 

2 

37.5 

*** 

82.0 

12.0 

7 

A 

189.5 

15030. 

2 

25.9 

*** 

77.1 

17.0 

8 

A 

189.5 

22620. 

2 

14.0 

**  * 

73.9 

21.0 

9 

A 

204.5 

440. 

2 

86.6 

*** 

95.3 

.6 

10 

A 

204.5 

825. 

2 

75.0 

*** 

91.5 

1.6 

11 

A 

204.5 

1200. 

2 

66.0 

*** 

88.8 

3.4 

12 

A 

204.5 

1500. 

2 

59.4 

*** 

86.4 

5.1 

13 

A 

204.5 

2040. 

2 

48.9 

*** 

83.0 

8.3 

14 

A 

204.5 

3060. 

2 

32.8 

*** 

77.8 

13.8 

15 

A 

204.5 

6060. 

2 

11.5 

*  *  * 

70.4 

22.5 

16 

A 

189.5 

36420. 

2 

4.5 

7.4 

70.5 

25.7 

17 

A 

204.5 

16020. 

2 

- 

3.1 

66.2 

28.6 

18 

A 

225.0 

3000. 

1 

- 

3.0 

66.0 

28.0 

19* 

A 

225.0 

3000. 

1 

- 

4.0 

66.0 

28.0 

20 

A 

245.0 

630. 

1 

— • 

4.0 

65.0 

27.0 

21* 

A 

245.0 

630. 

1 

- 

5.0 

65.0 

27.0 

22 

A 

265.0 

120. 

1 

- 

7.0 

65.0 

23.0 

23* 

A 

265.0 

120. 

1 

7.0 

65.0 

24.0 

24 

A 

285.0 

30. 

1 

~ 

11.0 

66.0 

19.0 

25* 

A 

285.0 

30. 

1 

- 

9.0 

66.0 

19.0 

26 

D 

189.5 

1020. 

1 

- 

- 

- 

80.0 

27 

D 

189.5 

3990. 

1 

- 

- 

- 

87.3 

28* 

D 

189.5 

3990. 

1 

- 

- 

- 

87.3 

29 

D 

189.5 

6780. 

1 

- 

- 

- 

87.5 

30 

D 

189.5 

8220. 

1 

- 

- 

- 

86.5 

31 

D 

189.5 

13260. 

1 

- 

- 

- 

88.5 

32 

D 

189.5 

14760. 

1 

- 

- 

- 

89.8 

33 

D 

204.5 

3480. 

1 

- 

- 

- 

87.5 

34 

D 

204.5 

5700. 

1 

- 

- 

- 

86.8 

35 

E 

189.5 

8880. 

1 

- 

- 

- 

91.9 

36* 

E 

189.5 

8880. 

1 

- 

- 

- 

92.0 

37 

E 

189.5 

14340. 

1 

- 

- 

- 

89.8 

38 

E 

189.5 

23400. 

1 

- 

- 

- 

89.7 

39* 

E 

189.5 

23400. 

1 

- 

- 

- 

88.5 

40 

E 

204.5 

5700. 

1 

- 

- 

- 

88.4 

41 

E 

204.5 

8100. 

1 

- 

- 

- 

87.9 

* 

Replicate  of  the  preceding  test. 

** 

is  the  number  of  independent  tests  combined  to  obtain  each  observation  y  . 

*** 

Originally  reported  but  not  observed;  see  text. 

-  No  value  reported. 
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Table  5.  Parameters  for  a-Pincne  Conversion 


Parameter 

Estimates 

4-Reaction 

for 

Model 

* 

Estimates 
5- Reaction 

for 

Model 

61 

-8.331 

.024 

-8.333 

4 

.025 

62 

-8.898 

+ 

.029 

-8.961 

+ 

.054 

03 

-8.242 

+ 

.341 

-8.196 

± 

.325 

04 

-5.389 

+ 

.081 

-5.438 

+ 

.087 

®5 

-11.945 

4- 

.698 

06 

19814. 

+ 

428. 

19785. 

+ 

457. 

67 

20828. 

+ 

474. 

20890. 

+ 

536. 

08 

17336. 

+ 

4079. 

17212. 

+ 

4203. 

09 

10321. 

+ 

915. 

10322. 

± 

918. 

0, 

19957. 

** 

10 

011 

269. 

+ 

83. 

279. 

± 

83. 

612 

-1976. 

+ 

64. 

-1985. 

± 

63. 

613 

-336. 

+ 

950. 

-259. 

± 

958. 

014 

-3873. 

± 

1624. 

-3781. 

± 

1555. 

°11 

.696 

4- 

.419 

.784 

± 

.492 

°21 

.000 

** 

.000 

Hit 

°22 

.391 

± 

.359 

.376 

± 

.348 

°31 

.358 

+ 

.412 

.426 

± 

.456 

°32 

.000 

** 

.000 

** 

°33 

.706 

± 

.426 

.732 

± 

.444 

°4i 

-.248 

± 

.344 

-.294 

.354 

°42 

.000 

** 

.000 

** 

°43 

-.504 

± 

.317 

-.493 

+ 

.314 

°44 

.744 

± 

.304 

.654 

+ 

.282 

95%  highest  posterior  density  intervals  calculated  from  Equation  (14). 

*• 

Posterior  estimates  were  not  obtained  for  these  parameters. 


l 


Expt • g 

u 

4- 

Table  6 

.  Final 

Residuals 

e  .  (0) 

Ul  ~ 

for  Example  2 . 

•Reaction 

Model 

5-Reaction  Model 

£  i 

ul 

(A) 

£u2 

(A+B) 

£u3 

(a+b+c) 

eu4 

<E) 

£ul 

(A) 

£u2 

(A+B) 

£u3 

(A+B+C) 

£u4 

(E) 

i 

-1.32 

_ 

-.37 

2.00 

-1.22 

- 

-.26 

1.69 

2 

-1.42 

- 

-.87 

1.10 

-1.32 

- 

-.76 

.79 

3 

.26 

- 

.24 

.88 

.43 

- 

.43 

.45 

4 

.28 

- 

-.15 

1.10 

.45 

- 

.06 

.72 

5 

.38 

- 

-.04 

.22 

.48 

- 

.12 

-.04 

6 

-1.13 

- 

.70 

-.81 

-1.12 

- 

.78 

-.96 

7 

-.32 

- 

-.26 

-.17 

-.43 

- 

-.29 

-.18 

8 

.66 

- 

.89 

-1.06 

.47 

- 

.74 

-.92 

9 

.88 

- 

.21 

.30 

1.00 

- 

.35 

-.11 

10 

.10 

- 

-.17 

.14 

.24 

- 

.04 

-.38 

11 

.31 

- 

-.07 

.16 

.45 

- 

.15 

-.34 

12 

.27 

- 

-.51 

.23 

.38 

- 

-.29 

-.23 

13 

-.04 

- 

-.86 

.42 

.01 

- 

-.69 

.08 

14 

-1.44 

- 

-1.56 

.90 

-1.52 

- 

-1.49 

.75 

15 

-.47 

- 

-1.47 

.63 

-.70 

- 

-1.61 

.77 

16 

.60 

.78 

.98 

-.36 

.44 

.72 

.74 

-.14 

17 

- 

-.12 

-.67 

.34 

- 

-.07 

-.87 

.50 

18 

- 

-.81 

.51 

-.48 

- 

-.76 

.38 

-.39 

19 

- 

.19 

.51 

-.48 

- 

.24 

.38 

-.39 

20 

- 

-.89 

.29 

-.56 

- 

-.88 

.22 

-.47 

21 

- 

.11 

.29 

-.56 

- 

.13 

.22 

-.47 

22 

- 

-.54 

-.31 

-.37 

- 

-.58 

-.32 

-.28 

23 

- 

-.54 

-.31 

.63 

- 

-.58 

-.32 

.72 

24 

- 

1.54 

.49 

-.15 

- 

1.51 

.58 

-.20 

25 

- 

-.46 

.49 

-.15 

- 

-.49 

.58 

-.20 

26 

- 

- 

- 

1.12 

- 

- 

- 

1.95 

27 

- 

- 

- 

-.92 

- 

- 

- 

-.61 

28 

- 

- 

- 

-.92 

- 

- 

- 

-.61 

29 

- 

- 

- 

-1.31 

- 

- 

- 

-1.16 

30 

- 

- 

- 

-2.37 

- 

- 

- 

-2.27 

31 

- 

- 

- 

-.41 

- 

- 

- 

-.42 

32 

- 

- 

- 

.90 

- 

- 

- 

.86 

33 

- 

- 

- 

.67 

- 

- 

- 

.72 

34 

- 

- 

- 

-.31 

- 

- 

- 

-.40 

35 

- 

- 

- 

1.26 

- 

- 

- 

1.17 

36 

- 

- 

- 

1.36 

- 

- 

- 

1.27 

37 

- 

- 

- 

.24 

- 

- 

- 

.16 

38 

- 

- 

- 

.91 

- 

- 

- 

.80 

39 

- 

- 

- 

-.29 

- 

- 

- 

-.40 

40 

- 

- 

- 

.42 

- 

- 

- 

.27 

41 

- 

- 

- 

.51 

- 

- 

- 

.35 

19 


All  i p  values  such  that  S(^)  <  65.50  lie  within  this  estimated  OS's  jo::.; 
confidence  region.  By  this  criterion,  the  model  with  k^  =  0  is  accei ‘ abl 
However,  as  indicated  earlier.  Equation  (14)  is  more  reliable  than  (13).  F 
this,  and  a  study  of  the  residuals,  we  conclude  that  the  5-reaction  model  i 
to  be  preferred. 
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APPENDIX:  DERIVATIVES  OF  S 


The  matrices  a  are  real  and  symmetric;  furthermore,  S  is  defined  only 
when  these  matrices  are  positive  definite.  The  following  derivative  relations 
then  hold: 
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The  relations  for  second  derivatives  follow  by  combination  of  (Al)  and  (A2) : 
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As  indicated  earlier,  if  response  h  is  absent  from  experiment  u  ,  the 
elements  °uhj  and  °ujh  are  rePlacec^  by  tlle  constant  dummy  values  5^  . 
Note  also  that  the  symmetry  of  has  been  used  to  express  these  deriva¬ 

tives  in  terms  of  elements  on  and  below  the  diagonal. 

The  derivatives  required  for  Equation  (9)  are  obtained  as  follows: 
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Equations  (A6) ,  (hi),  and  (a9)  evaluated  at  0^  and  o  provide  the  coef¬ 
ficient  matrix  A  of  the  normal  equations.  Equations  (A5)  and  (A8)  give  the 
right-hand  column  vector. 


The  residuals  e  .  and  c  .  are  expressed  as  functions  of  0  by  use 
ur  ug  , 


of  Equation  (1).  The  0-derivative  in  Equation  (A6)  is  expanded  to  give: 
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The  second-derivative  terms  are  unimportant  if  the  data  are  well  fitted; 

compare  Solutions  1  and  2  in  Table  3. 

If  the  experiments  have  different  weights  wy  as  in  Table  4,  then 

e  .  e  .  and  its  derivatives  should  be  multiplied  by  w  throughout  the 
ui  uj  J  u 

development.  As  usual,  the  matrix  a  is  defined  for  experiments  of  unit 
weight. 
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